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Abstract 

Subject of this paper is an implementation of a well-known Motzkin- 
Burger algorithm, which solves the problem of finding the full set of solu- 
tions of a system of linear homogeneous inequalities. There exist a number 
of implementations of this algorithm, but there was no one in Maple, to 
the best of the author's knowledge. 

1 The problem 

Consider a linear homogeneous system of inequalities with rational coefficients 

lj{x) = GjiXi + . . . + ajnXn S^Q, ttj^ £ Q, (j = l,...,m). (1) 

Let r be the rank of the matrix of 

Recall some common facts from the theory of convex sets (see for example, 
^). Every solution set Xi, . . . , Xp of a system generates a convex cone C 
of solutions 

kiXi + ... + kpXp, h ^ 0. (2) 

We call a finite set of solutions Xi, . . . , Xp e Q" a set of generators for C 
if every element of C is a conical combination ^ oi Xi,. . . ,Xp. Hereinafter 
we will consider the minimal sets of generators only, i.e. such that none of the 
Xi,i — 1,. . . ,p is expressed from the others with positive coefficients. We call 
d the dimension of a cone C if d is the dimension of the minimal linear space 
spanned by C. 

A solution locus of each inequality is a half-space of Q". For any subsystem 
{lj{x) ^ 0,j € J} of its solution set is a convex polyhedral cone. The 
faces of the latter are the intersections of C with solution loci of some equation 
subsystems {lj(x) = 0, j G / C J} with linearly independent left-hand sides. 

The problem discussed in this paper is to solve Q a over Q". This problem 
is closely related to the following one. Extend the system |^ by a new inequality 

l{x) = aiXi -I- ... -I- ttnXn ^ 0. (3) 

Not all the points of the cone C satisfy lj2Jl, thus the system ^ together with 
^ define a new cone C*. We are interested in computing the set of generators 
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for C* given the set of generators for C. The Motzkin-Burger algorithm solves 
this problem. 

The iteration of the Motzkin-Burger algorithm solves the first problem by 
means of extending the trivial system {0 ^ 0} by the inequalities of Q one 
by one. 

2 Algorithm description 
2.1 Theoretical study 

Let us recall some well-known facts (generally following 0, 0). 

The set of generators for a cone C* consists of generators for cone C (which 
satisfy the inequality l{x) < 0) and the whole set of generators for cone C H H 
where H := {l{x) = 0}. Thus, in order to compute the cone C* we need to 
study both the structure of the cone C and structure of the intersection. 

The cone C is a Minkovsky sum of its linear subspace L and the strongly 
convex cone P (that is the cone with apex at the origin and contain no line 
through the origin). Vectors u of the base U oi L satisfy {lj(u) = 0,j = 
1, . . . , m}, whereas any element v of set V of generators of P satisfies lj{v) < 
for at least one j. 

Let us denote vectors X belonging to C/ or y by X+ if 1{X) > 0, if 
1{X) < and if 1{X) — 0. Also denote the linear subspace of the new 
cone with its base and the strongly convex cone of the new cone with its set of 
generators by L* , U*, P*, V* , respectively. Hence the vector u*^ belongs to U*, 
while u~ , v'^ and belong to V* . Although, {u E U : l{u) = 0} is not U*. By 
discarding all it+-s and v~^-s we lose a number of generators for the cone C and 
have to be replaced them by the whole set of generators for H CiC. 

Since H DC ~ {H n i) U [H n P) we can describe how to convert U to U* 
and V to V* separately. 

There are two cases with respect to the intersection HDL: H D L — L oi 
Hr\L^ L. In the first case, U* = U. 

In the second case the new inequality reduces the dimension of L by one. 
Choose one u~ € U (if there is no such vector in U, take — Let us form 
diniL — 1 pairs of vectors and Uj G U \ {u^}. For every pair consider the 
conic combinations 

au^ + buj, a, 6 ^ 0. (4) 
The intersection of H with the cone Q is the ray of positive multiplies of 

u* :— —l(yr)uj + l{uj)u^ . (5) 

It is clear that l{u*) = 0. Thus, dimi — 1 vectors u* are the base of L*. 
Note that the representation of the base of L is not unique since we can choose 
different vectors for u~ if they exist. 

The conversion of V to V* also depends on existence of abovementioned u~ . 
The element u~ belongs to V* if it exists. New inequality lO induces the affine 
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transformation 

-l{u^)v + l{v)u^ (6) 

of elements v ^ V that results in that every vector © satisfies 0. It is 
evident that the set {u^ , —l{u^)v + l{v)u~ for all v G V} is minimal in the 
abovementioned sense, thus the latter is the set of generators V* . 

If there exist no such u~ £ C/, the transformation lO of is impossible. 
The new inequality separates P into two parts. As already was mentioned, the 
elements v~ belong to V* . 

To find all the vectors that lie in H let us consider any pair w^, t)+ e V . For 
this pair the combination 

vl:= -l{v~)v+ +l{v+)v~ (7) 

lies in H. The set of all u^-s generates the convex polyhedral cone in H , but 
there are too many elements of {w^} for it to be the set of generators for the 
latter cone. The minimal faces that contain a pair of generators are exactly 
two-dimensional ones. Therefore, in order to avoid the superfluous solutions 
we only need to reject all the pairs of generators that are not lying on faces of 
dimension 2. 

Checking if two vectors lie on a face of dimension 2 may be done in two ways. 
The first way is to check whether there exists {lj{x),i G J,\J\ = r — 2} with 
linear independent lj{x) such that lj{v~) — lj{v^) — holds true for all j G J. 
The second way is to check there exist no third v E V such that lj{v) = for 
all j such that lj{v^) ~ lj{v^) = (are not taken into account the rank r nor 
linear independence of Ij-s). 

We therefore come to the algorithm that solves a system of linear inequalities 
over Q": 



Input: S := {li{x), . . . , lm{x)} — the left-hand sides of the inequalities 

of J^l, n — the space dimension. 
Output: U — {ui, . . . , ut} — the base of the maximal linear subspace L of 

the solution cone of 

V = {vi, . . . ,Vs} — the set of generators for the strongly convex 
cone P of the solutions of the system Q . 



{(1,...,0),...,(0,. ..,!)}; 

0; 

{0 0}; 



1 • U current 

2- ^current ■ 

3- Scurrent 

4. i 1; 

5. I :— li{x); 

6. if Elu G Ucurrent- l{u) ^ 0, then 

^current • — ^(Wj)'U, Ui G Ucurrent} ^ 

n :— —{l{u)/\l{u)\)u; 

Vcurrent := {?^, -l{n)Vi + l{Vi)n, Vi G Vcurrent); 



7. else 



yLrrent - ^ ^4«rre„t | l{v) ^ 0}; T^l, 
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for y{vk,Vs) G Klrrent : K^k) < 0, > do 

S* := {lj{x) e Scurrentl l]{vk) = l]{Vs) = 0}; 

if 5 V then 

for V W e Vcurrent \ {vk, Vs} 

for V/j (a;) £ S* do 
if lj{v) ^ then 

^current — Klrrent ^ + /(wfc)Ws}; 

end if; 
end do; 
end do; 
end if; 
end do; 
end if; 

8- Vcurrent '■= Vcurrent ^ ^current'i 

9- ^current • — ^current U {^«(*^)}; 

10. S:=S\{kix)}; 

11. + 

12. if S* = then 

goto 14; 
end if; 

13. goto 5; 

14. return Ucurrent^ Vcurrent • 



2.2 Some improvements 

There are some tricks concerning the preparation of the system in order to lower 
the number of inequahties and/or variables. 

Firstly, the system may have no occurrences of some of Xfe, fc = 1, . . . , n. In 
this case we can "clean up" the system and thus reduce the number of variables. 

Secondly, we can avoid the linear subspace L of cone of solutions by perform- 
ing the change of variables. This trick also enables one to simplify the original 
system of inequalities. 

The dimension of i is n — r, where r is the rank of matrix of as was 
mentioned above. Let us define new r variables j/i, . . . , j/^ to equal any r linear 
independent left-hand sides of (for example, first r ones): 

lj{x)^~y, i^O, j = l,...,r, (8) 
lj{x)s^O, j = r + l,...,m. (9) 

By solving the system {lj{x) — ~yj,j = 1, . . . , r} for {xi, . . . , Xn} we obtain 
the (ri — r)-dimensional space of solutions 
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(10) 

■^ir — fr (j/l : ■ ■ ■ :yr: ^ir+i 7 ■ ■ ■ 7 •'^^,^)^ 

(11) 

where (ii, . . . , in) is the permutation of indices 1, . . . , rt (we don't know a priori 
which variables most hkely to be solved for). 

Upon substitution I|1Q|I the system does not depend on Xi^_^^ , ■ ■ ■ , Xi^ be- 
cause /r+i(a;)7 • • ■ , lm{x) are the linear combinations of the first r ones. There- 
fore, substitution pU|) reduces the system to that of m inequalities for r indepen- 
dent variables t/i, . . . , y^- The r inequalities of reduced system are simplified to 
be just —iji ^ 0. This inequality subsystem has the a priori known solution E^. 
— the r-dimensional linear space base. Hence we may efface these inequalities 
from the reduced system. More precisely, solving the system we iterate Motzkin- 
Burger algorithm starting from the system of inequalities {—yi ^ 0, « = 1, . . . , r} 
and U = % and V = Er- Thus the system is reduced to of m — r inequalities for 
r independent variables. Because ?7 = at the start of algorithm, there is no 
need in any parts of algorithm except for the most complicated part that is the 
conversion of V to V* when there is no u" exists. 

The application of the Motzkin-Burger algorithm yields a number of so- 
lutions of the reduced system. By means of n — r substitutions (|10|l and 
(xi^^j, . . . , Xi^) — E^^j., where are the base vectors of (n — r)-dimensional 

linear space, we then obtain the solution set of the original system. 

The case of r = n is also of interest for simplifying the system despite that 
there is no lowering the number of variables. The system is reduced to that 
of m — n inequalities for n independent variables in this case. The calculation 
procedure is repeated in full except for that there is no need to consider the 
substitution the {n — r)-dimensional linear space base vectors to the several 
components {xi^_^-^ , ■ ■ ■ , Xi^) of any solution. 

Also of interest is the case of r = to where m is the number of inequalities. 
Upon the abovementioned substitution, the system is reduced to the diagonal 
one. Therefore, the whole system have a priori known solution in terms of new 
variables. 

We therefore come to the algorithm of solving the system ^ with special 
preparation prior to the application of the Motzkin-Burger algorithm: 

Input: S := {/i(a;), . . . , lm.{x)} — the left-hand sides of the inequalities 

of 1^, n — the space dimension. 
Output: U = {ui, . . . , ut} — the base of the maximal linear subspace L of 

the solution cone of 
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V = {vi, . . . , Vs} — the set of generators for the strongly convex 
cone P of the solutions of the system 

1. Find all indices Bad C {I, . . . ,n} such that S has no occurrences of Xj, 
j G Bad; 

2. Renumber the indices (not encountered in Bad) of variables so that S 
explicitly depend on n — \Bad\ variables only; 

3. Choose Base — the maximal linear independent subsystem of {lj{x),j = 
1, . . . , to}; r := |_Base|; 

4. Solve {Basci — —yi, i — 1, . . . , r} for {xi, . . . , Xn} to obtain set of identi- 
ties (fTU)) : 

5. / := {v+i, . . . , «„}; 

6. U' := 0; 

7. V :={(1,0,...,0),...,(0,...,0,1)}; 

" V ' 

r components 

8. if r < m then 

(a) Substitute ifTUI) to S in order to obtain S"; 

(b) Reorder the inequalities S such that —yi ^ are the first r ones. 

(c) Iterate the part of Motzkin-Burger algorithm mentioned above ex- 
cluding the items 1, 2, 3, 4, 6, considering Ucurrent = U* , Vcurrent = 
V', S current = S' , i — r + 1 as the initial condition. 

end if; 

9. E:={ (1,0, . .., 0) ,...,(0,...,0,1)}; 

n — r components 

10. For each element e G _B substitute {yi — ... = yr = 0} and (x/^, . . . , 
a;/„_r) — e into (|10|l and form u :— (xi, . . . , G U ; 

11. For each element v G V substitute {(?/i, . . . , yr) — v} and (x/j = . . . = 
^in-r) — iiito H1Q(I and form v := {xi, . . . ,a;„) G V; 

12. return U, V. 

2.3 Algorithm complexity 

Let us not distinguish arithmetic and comparison operations. Let n be the 
number of independent variables, p be the number of elements in V at the 
current step, I be the number of , v £ V , fc be the number of ,v E V, q be 
the number of inequalities already examined by the current step and to be the 
number of inequalities in total. 

The calculation of the value of l{x) has the complexity 2n — 1. If there exists 
such u € U that l{u) 0, than the conversion U U* and V ^ V* requires 
the calculation of n — 1 differences of the kind oil{u^)uj — l{uj)u^ with the total 
complexity 2 • 2n + 1 = 4n + 1. Therefore, the overall complexity of converting 
both lists is 2(n - l)(4n -h 1) ~ 8n^. 

The case when no such u G U exists is of baffling complexity. Let us estimate 
the complexity of one iteration of this case as a function f{p, q, n). For efficiency 
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reasons wc suppose that all the values of Ij (x) G S for all v G V are computed 
beforehand in order to exclude repeated calculations. 

Selection of v~ takes I operations. Next step — to select the pairs v~,v~^ 
— requires 2kl more operations. The number of pairs reaches the maximum if 

2 

A; = / = I SO let fc/ = ^ below. 

For every pair v~,v~^ we need to choose those inequalities lj{x) for which 
lj{v~) = lj{v^) is true. This requires 2q operations. Then we should examine 
every pair on whether p — 2 elements of V do not zero all the chosen inequal- 
ities. Let the number of the chosen inequalities be as large as possible i.e., |. 

2 3 2 

Therefore, checking all the pairs have the complexity ^(p — 2)| = ^ — 

Let all the pairs satisfy the conditions above, then the overall complexity for 

all new the elements v is 

In total, f{p, q, n) = pq{2n - 1) + R + + - + ^ = Ipip^q + lOp - 

2pq + 16gn - 8g + 4) ~ \p^q. 

These intensive calculations are possible since po = 3 because up to p 2 

there is only one (or no one) pair of generators for V. Therefore, up to one new 

generator is produced by algorithm for it to replace another. 

2 

As was mentioned, number p^ of elements v £V grows as + ^-^f^ on ev- 
ery fc-th step. Iterating this function one can see that pfe = Pk{Pk-i) = 0{p'q ^ ), 
so the worst-case complexity of pure iteration of the most complicated part of 
algorithm is f{po,m,n) = 0{mp^ ^ /8) = 0(mpg^'"' = O(rnpo^"') = 
0(m32""). 

Practice shows that, in fact, such a terrible complexity is almost impossible 
to happen. There are many pairs rejected on every step. The number of pairs 
generally is not too large. Number of probe inequatilities for every pair usually 
lower than |. Usually systems of full rank are incompatible for sufficiently large 
number of inequalities. Overall computation time nevertheless depends on n for 
the moderate systems. Despite of it, the data grows exponentially. 

2.4 Practical experience 

We implemented abovementioned algorithms in Maple as one package 
motzkin_burger consists of three user procedures: Conehull, MB and 
CheckSolutions. The first of these solves the problem to compute a set of 
generators for solution cone of system of inequalities. The second is the single 
Motzkin-Burger iteration. Last of these is the tool that allow user to check up 
the set of solutions on a correctness. Prototypes of these are 

Conehull (L , x , n , options) 

MB(L,U,V,x,n) 

CheckSolutions (L , S ,x , n) 
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Here L is a list of homogeneous polynomials of degree 1 in terms of variables 
Xi, n is dimension of space of solutions; U is the base of linear space of cone 
of solutions, V is the set of generators of strongly convex cone in latter cone; S 
is any solution set. Parameter options, at this moment, may accept only one 
value, "as is" , what gives an instruction do not perform change of variables. 

Both Conehull and MB return an exprseq that consists of two 2-dimensional 
lists which are abovementioned lists U and V. CheckSolutions returns an 
exprseq too, but it might by interpret in another way: first is the list of vectors 
that are true solutions whereas the second is the list of uncorrect solutions. 

The timings presented in the table below are obtained in Maple 7 on com- 
puter based on Duron 700Mhz, 256 RAM. For every combination of (n,m,r) 
there was computed random system having these parameters, with sufficiently 
small integer coefficients. Here n is space dimension, m is the number of in- 
equalities in the system and r is the rank of the system. Values ti and t2 are 
times in seconds for solving the latter systems using or not optional parame- 
ter of Conehull procedure. All these systems are saved in attached to paper 
package codes text files. 



space dimen- 


number of in- 


rank, r 


Conehull, option 


Conehull, 


sion, n 


equalities, m 




"as is", ti, sec 


t2, sec 


5 


5 


5 


0.010 


0.120 


5 


7 


3 


0.090 


0.050 


10 


10 


10 


0.130 


0.261 


10 


15 


5 


0.150 


0.180 


20 


20 


20 


1.783 


2.103 


20 


30 


10 


1.512 


1.021 


30 


30 


15 


2.864 


2.003 


40 


40 


20 


8.663 


4.316 


40 


40 


30 


45.326 


17.775 


50 


50 


40 


89.118 


1816.392 


50 


50 


45 


73.766 


1549.628 



One can sec that Conehull with option "as is" computes this examples 
with always increasing time (with rare exceptions). Without the option it be- 
have more complicated: for every fixed n and m the case of the systems of full 
rank or "almost" full rank may be computed more faster than in the case of 
small rank. Of course, the option slows down the performance (this is the con- 
sequence of inefficient Maple subs implementation we are using) starting from 
some dimension, but we hope that the systems of larger dimension will be com- 
puted more faster with this option. Nevertheless, to improve this performance 
bottleneck is the main aim of further work. 
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